# Replication File for "Corruption Victimization and Anti-Incumbent Voting"
# Aimee Bourassa, Rebecca Weitz-Shapiro, and Matthew S. Winters
# Governance

# This file runs the analyses reported in "Corruption Victimization and 
 # Anti-Incumbent Voting."  The replication data is based on cleaned 
 # and recoded versions of the 2014 AmericasBarometer data (as included in
 # the 2004 - 2015 Grand Merge File found on the LAPOP website); the 
 # 2016 (Round 6) Afrobarometer data; and country-level data from the Polity
 # project, the World Development Indicators, and the QoG Expert Survey II.

# Variables in replication datasets:
 # ccode_abc - three-letter country identifier
 # vote_incumbent - 0/1 indicator for intention to vote for incumbent in next national election
 # vote_yes - 0/1 indicator for intention to turnout to vote in next national election
 # bribe_victim - 0/1 indicator for exposure to corruption 
 # female - 0/1 indicator for female
 # age - age in years
 # edu - 0...6 categorical variable: (0) no formal education; (1) incomplete primary school;
  # (2) complete primary; (3) incomplete secondary; (4) complete secondary; (5) some university
  # or technical/post-secondary education; (6) completed university education and beyond
 # quantSES - respondent's in-country quintile on a household asset index (AmericasBarometer)
  # or the Lived Poverty Index (Afrobarometer)
 # party - 0/1 indicator for partisanship with the incumbent party
 # wt - AmericasBarometer survey weight
 # estratopri - AmericasBarometer PSU indicator
 # withinwt - Afrobarometer survey weight
 # p_polity2 - Polity score
 # gdppcgrowth - GDP per capita growth rate from year of survey from WDI
 # inflation - inflation rate from the year of the survey from WDI
 # politicized - measure of bureaucratic politiciziation from QoG Expert Survey II;
  # seven-point scale ranging from "hardly ever" (1) to "almost always" (7)
 # closed4 - index based on four questions from QoG Expert Survey II asking about 
  # the extent to which "(1) "public sector employees are hired via a formal 
  # examination system"; (2) "once one is recruited as a public sector employee, 
  # one remains a public sector employee for the rest of one's career"; 
  # (3) "the terms of employment for public sector employees are regulated by 
  # special laws that do not apply to private sector employees"; and
  # (4) "senior public officials are recruited from within the ranks of the 
  # public sector"; each question is originally scored on a seven-point scale
  # from "hardly ever" (1) to "almost always (7); we take the average
 # efficient - measure of the extent to which "public sector employees strive 
  # to be efficient" from the QoG Expert Survey II; seven-point scale ranging
  # from "hardly ever" (1) to "almost always" (7)
 # absent - measure of the extent to which "public employees are absent from 
  # work without permission" from the QoG Expert Survey II; seven-point scale
  # ranging from "hardly ever" (1) to "almost always" (7)

##############
# BACKGROUND #
##############

library(readstata13)
library(stargazer)
library(estimatr)
library(survey)
library(dplyr)
library(ggplot2)
library(lme4)
library(interplot)
library(gridExtra)
library(R2HTML)

# centering function
center_vector <- function(x){
  as.vector(scale(x, scale=F))
}

########################
# LOAD AND FORMAT DATA #
########################

lapop <- read.dta13("lapop_replication.dta")
afro <- read.dta13("ab_replication.dta")
country <- read.dta13("countrylevel_replication.dta")

# Log Transformation
country$loggdppcppp <- log(country$gdppcppp)

vars <- c("vote_incumbent", "vote_yes", "bribe_victim", "female", "age", "edu", 
  "quantSES", "ccode_abc", "party", "wt", "some_contact", "bribe_police", 
  "bribe_school", "bribe_health", "bribe_courts", "bribe_doc")
small.lapop <-lapop[,vars]
small.afro <- afro[,vars]

# Create a Country-Level Corruption Prevalence Measure
lapop.bribe.avg <- aggregate(lapop$bribe_victim ~ lapop$ccode_abc, FUN="mean")
colnames(lapop.bribe.avg) <- c("ccode_abc", "bribe_victim_avg")
afro.bribe.avg <- aggregate(afro$bribe_victim ~ afro$ccode_abc, FUN="mean")
colnames(afro.bribe.avg) <- c("ccode_abc", "bribe_victim_avg")

bribe.avg <- rbind(lapop.bribe.avg, afro.bribe.avg)

country <- merge(country, bribe.avg, by="ccode_abc")

# Merge in Country-Level Data
small.lapop2 <- merge(small.lapop, country, by="ccode_abc")
small.afro2 <- merge(small.afro, country, by="ccode_abc")
small.lapop2$survey <- "lapop"
small.afro2$survey <- "afro"

combined.data.mlm <- as.data.frame(rbind(small.lapop2, small.afro2))

# Centering Variables 
 # Note: needs to be done for each dataset so that the centering is within-sample

# politicization
small.afro2$politicized_zero <- center_vector(small.afro2$politicized)
small.lapop2$politicized_zero <- center_vector(small.lapop2$politicized)
combined.data.mlm$politicized_zero <- center_vector(combined.data.mlm$politicized)
# closed4
small.afro2$closed4_zero <- center_vector(small.afro2$closed4)
small.lapop2$closed4_zero <- center_vector(small.lapop2$closed4)
combined.data.mlm$closed4_zero <- center_vector(combined.data.mlm$closed4)
# efficiency
small.afro2$efficient_zero <- center_vector(small.afro2$efficient)
small.lapop2$efficient_zero <- center_vector(small.lapop2$efficient)
combined.data.mlm$efficient_zero <- center_vector(combined.data.mlm$efficient)
# absent
small.afro2$absent_zero <- center_vector(small.afro2$absent)
small.lapop2$absent_zero <- center_vector(small.lapop2$absent)
combined.data.mlm$absent_zero <- center_vector(combined.data.mlm$absent)
# polity 
small.afro2$p_polity2_zero <- center_vector(small.afro2$p_polity2)
small.lapop2$p_polity2_zero <- center_vector(small.lapop2$p_polity2)
combined.data.mlm$p_polity2_zero <- center_vector(combined.data.mlm$p_polity2)
# loggdppc
small.afro2$loggdppcppp_zero <- center_vector(small.afro2$loggdppcppp)
small.lapop2$loggdppcppp_zero <- center_vector(small.lapop2$loggdppcppp)
combined.data.mlm$loggdppcppp_zero <- center_vector(combined.data.mlm$loggdppcppp)
# cce
small.afro2$cce_zero <- center_vector(small.afro2$cce)
small.lapop2$cce_zero <- center_vector(small.lapop2$cce)
combined.data.mlm$cce_zero <- center_vector(combined.data.mlm$cce)
# bribe-level
small.afro2$bribe_victim_avg_zero <- center_vector(small.afro2$bribe_victim_avg)
small.lapop2$bribe_victim_avg_zero <- center_vector(small.lapop2$bribe_victim_avg)
combined.data.mlm$bribe_victim_avg_zero <- center_vector(combined.data.mlm$bribe_victim_avg)

###########
# TABLE 1 #
###########

predictors <- c("bribe_victim", "female", "age", "edu", "quantSES")
predictors.f <- paste(predictors, collapse= "+")

fit.combined.pooled <- lm(paste("vote_incumbent ~", predictors.f), data=combined.data.mlm)
fit.combined.fe <- lm(paste("vote_incumbent ~", predictors.f, "+as.factor(ccode_abc)"), 
  data=combined.data.mlm)

stargazer(fit.combined.pooled, fit.combined.fe,
  se = starprep(fit.combined.pooled, fit.combined.fe),
  column.labels=c("Pooled", "Fixed Effects"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES"),
  keep.stat = c("n", "rsq"), omit="ccode_abc", notes="Robust standard errors",
  type="html", out="pooledmodels.html")

########################
# LOGISTIC REGRESSIONS #
#    (FOOTNOTE 27)     #
########################

fit.combined.pooled.logit <- glm(paste("vote_incumbent ~", predictors.f), 
  data=combined.data.mlm, family="binomial")
fit.combined.fe.logit <- glm(paste("vote_incumbent ~", predictors.f, "+as.factor(ccode_abc)"), 
  data=combined.data.mlm, family="binomial")

stargazer(fit.combined.pooled.logit, fit.combined.fe.logit,
  column.labels=c("Pooled", "Fixed Effects"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES"),
  keep.stat = c("n", "rsq"), omit="ccode_abc", notes="Logistic regression models",
  type="html", out="pooledmodels_logit.html")

##############################
# TABLE 2: MULTILEVEL MODELS #
##############################

fit.mlm.pol <- lmer(vote_incumbent ~ bribe_victim + female + age + edu + quantSES + 
  (1+bribe_victim | ccode_abc) + bribe_victim*politicized_zero +
  bribe_victim*p_polity2_zero + bribe_victim*loggdppcppp_zero, 
  data= combined.data.mlm)                       
fit.mlm.closed <- lmer(vote_incumbent ~ bribe_victim + female + age + edu + quantSES + 
  (1+bribe_victim | ccode_abc) + bribe_victim*closed4_zero +
  bribe_victim*p_polity2_zero + bribe_victim*loggdppcppp_zero, 
  data= combined.data.mlm)                       
fit.mlm.eff <- lmer(vote_incumbent ~ bribe_victim + female + age + edu + quantSES + 
  (1+bribe_victim | ccode_abc) + bribe_victim*efficient_zero +
  bribe_victim*p_polity2_zero + bribe_victim*loggdppcppp_zero, 
  data= combined.data.mlm)                       
fit.mlm.absent <- lmer(vote_incumbent ~ bribe_victim + female + age + edu + quantSES + 
  (1+bribe_victim | ccode_abc) + bribe_victim*absent_zero +
  bribe_victim*p_polity2_zero + bribe_victim*loggdppcppp_zero, 
  data= combined.data.mlm)                       

stargazer(fit.mlm.pol, fit.mlm.closed, fit.mlm.eff, fit.mlm.absent,
  column.labels=c("Politicization", "Closedness", "Efficiency", "Absenteeism"),
  dep.var.labels.include=FALSE, dep.var.caption="", model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", "SES", 
   "Politicization","Closedness","Efficiency","Absenteeism",
   "Polity2", "Log(GDP per capita)", 
   "Bribe Victim*Politicization","Bribe Victim*Closedness",
   "Bribe Victim*Efficiency","Bribe Victim*Absenteeism",
   "Bribe Victim*Polity2", "Bribe Victim*Log(GDP per capita)"),
   keep.stat = c("n"), type="html", out="varyingslopemodels.html")

##########################################
# FIGURE 1: PLOTS FROM MULTILEVEL MODELS #
##########################################

pol.plot <- interplot(fit.mlm.pol, var1="bribe_victim", var2="politicized_zero") + theme_bw() +
  labs(x="Politicization", y="Marginal Effect of Bribe Victimization", title="") + ylim(-.2,.025) + geom_hline(yintercept = 0, linetype="dashed")
closed.plot <- interplot(fit.mlm.closed, var1="bribe_victim", var2="closed4_zero") + theme_bw() +
  labs(x="Closedness", y="", title="") + ylim(-.2,.025) + geom_hline(yintercept = 0, linetype="dashed")
eff.plot <- interplot(fit.mlm.eff, var1="bribe_victim", var2="efficient_zero") + theme_bw() +
  labs(x="Efficiency", y="Marginal Effect of Bribe Victimization", title="") + ylim(-.2,.025) + geom_hline(yintercept = 0, linetype="dashed")
abs.plot <- interplot(fit.mlm.absent, var1="bribe_victim", var2="absent_zero") + theme_bw() +
  labs(x="Absenteeism", y="", title="") + ylim(-.2,.025) + geom_hline(yintercept = 0, linetype="dashed")

allplots2 <- grid.arrange(pol.plot, closed.plot, eff.plot, abs.plot,
  ncol=2)
  
ggsave("mlmplots.png",allplots2, width=9, height=9)

################################
# TABLE 3: PARTISAN MODERATION #
################################

fit.party.combined.pooled <- lm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + party + bribe_victim*party, data=combined.data.mlm)
fit.party.combined.fe <- lm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + party + bribe_victim*party + as.factor(ccode_abc), data=combined.data.mlm)

stargazer(fit.party.combined.pooled, fit.party.combined.fe,
  se= starprep(fit.party.combined.pooled, fit.party.combined.fe),
  column.labels=c("Pooled", "Fixed Effects"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES", "Incumbent Copartisan", "Bribe*Copartisan"),
  keep.stat = c("n", "rsq"), omit="ccode_abc",
  type="html", out="pooledmodels_partisanship.html")

#########################################
# APPENDIX B: BRIBE VICTIMIZATION RATES #
#########################################

lapop.bribe.avg.2 <- lapop.bribe.avg
afro.bribe.avg.2 <- afro.bribe.avg
lapop.bribe.avg.2[,2] <- round(lapop.bribe.avg.2[,2], 2)
afro.bribe.avg.2[,2] <- round(afro.bribe.avg.2[,2], 2)
HTML(lapop.bribe.avg.2, "BribeRates_LAPOP.html")
HTML(afro.bribe.avg.2, "BribeRates_Afro.html")

################################## 
# APPENDIX C: COUNTRY-BY-COUNTRY #
##################################

# SEE BELOW

#########################
# APPENDIX D: BY SURVEY #
#########################

# APPENDIX TABLE 1

fit.lapop.pooled <- svyglm(paste("vote_incumbent ~", predictors.f), 
  design = svydesign(ids=~1, strata=lapop$estratopri, weights=lapop$wt, data=small.lapop))
fit.lapop.fe <- svyglm(paste("vote_incumbent ~", predictors.f, "+as.factor(ccode_abc)"), 
    design = svydesign(ids=~1, strata=lapop$estratopri, weights=lapop$wt, data=small.lapop))
fit.ab.pooled <- svyglm(paste("vote_incumbent ~", predictors.f), 
  design = svydesign(ids=~1, strata=NULL, weights=afro$wt, data=small.afro))
fit.ab.fe <- svyglm(paste("vote_incumbent ~", predictors.f, "+as.factor(ccode_abc)"), 
    design = svydesign(ids=~1, strata=NULL, weights=afro$wt, data=small.afro))

stargazer(fit.lapop.pooled, fit.lapop.fe, 
    fit.ab.pooled, fit.ab.fe,
  column.labels=c("LAPOP<br>Pooled", "LAPOP<br>FE", "Afrobarometer<br>Pooled", 
    "Afrobarometer<br>FE"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES"),
  notes = "OLS regressions estimated with survey weights.",
  keep.stat = c("n", "rsq"), omit="ccode_abc",
  type="html", out="pooledmodels_bysurvey.html")

#######################
# APPENDIX E: TURNOUT #
#######################

fit.t.combined.pooled <- lm(paste("vote_yes ~", predictors.f), data=combined.data.mlm)
fit.t.combined.fe <- lm(paste("vote_yes ~", predictors.f, "+as.factor(ccode_abc)"), 
  data=combined.data.mlm)
length(coef(fit.t.combined.fe)) - length(coef(fit.t.combined.pooled)) + 1 # calculate J
fit.t.lapop.pooled <- svyglm(paste("vote_yes ~", predictors.f), 
  design = svydesign(ids=~1, strata=lapop$estratopri, weights=lapop$wt, data=small.lapop))
fit.t.lapop.fe <- svyglm(paste("vote_yes ~", predictors.f, "+as.factor(ccode_abc)"), 
    design = svydesign(ids=~1, strata=lapop$estratopri, weights=lapop$wt, data=small.lapop))
length(coef(fit.t.lapop.fe)) - length(coef(fit.t.lapop.pooled)) + 1 # calculate J
fit.t.ab.pooled <- svyglm(paste("vote_yes ~", predictors.f), 
  design = svydesign(ids=~1, strata=NULL, weights=afro$wt, data=small.afro))
fit.t.ab.fe <- svyglm(paste("vote_yes ~", predictors.f, "+as.factor(ccode_abc)"), 
    design = svydesign(ids=~1, strata=NULL, weights=afro$wt, data=small.afro))
length(coef(fit.t.ab.fe)) - length(coef(fit.t.ab.pooled)) + 1 # calculate J

# APPENDIX TABLE 2: Turnout

stargazer(fit.t.combined.pooled, fit.t.combined.fe, fit.t.lapop.pooled, 
  fit.t.lapop.fe, fit.t.ab.pooled, fit.t.ab.fe,
  column.labels=c("All Pooled", "All FE", "LAPOP<br>Pooled", "LAPOP<br>FE", 
    "Afrobarometer<br>Pooled", "Afrobarometer<br>FE"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES"),
  notes = "OLS regressions estimated with survey weights.",
  keep.stat = c("n", "rsq"), omit="ccode_abc",
  type="html", out="turnout_pooledmodels_bysurvey.html")

##################################################
# APPENDIX F: ACROSS DEMOCRACIES/NON-DEMOCRACIES #
##################################################

fit.combined.pooled.dem <- lm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES, data=subset(combined.data.mlm,p_polity2>=6))
fit.combined.fe.dem <- lm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + as.factor(ccode_abc), data=subset(combined.data.mlm,p_polity2>=6))
fit.combined.pooled.nondem <- lm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES, data=subset(combined.data.mlm,p_polity2<6))
fit.combined.fe.nondem <- lm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + as.factor(ccode_abc), data=subset(combined.data.mlm,p_polity2<6))

# APPENDIX TABLE 3

stargazer(fit.combined.pooled.dem, fit.combined.fe.dem,
 fit.combined.pooled.nondem, fit.combined.fe.nondem,
  se = starprep(fit.combined.pooled.dem, fit.combined.fe.dem,
 fit.combined.pooled.nondem, fit.combined.fe.nondem),
  column.labels=c("All Data<br>Pooled<br>Dem", "All Data<br>FE<br>Dem",
    "All Data<br>Pooled<br>Non-Dem","All Data<br>FE<br>Non-Dem"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES"),
  keep.stat = c("n", "rsq"), omit="ccode_abc", note="Robust standard errors.",
  type="html", out="pooledmodels_dem_non.html")

######################################
# APPENDIX G: PARTISANSHIP BY SURVEY #
######################################

fit.party.lapop.pooled <- svyglm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + party + bribe_victim*party, design = svydesign(ids=~1, strata=lapop$estratopri, 
    weights=lapop$wt, data=lapop))
fit.party.lapop.fe <- svyglm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + party + bribe_victim*party + as.factor(ccode_abc), 
    design = svydesign(ids=~1, strata=lapop$estratopri, 
    weights=lapop$wt, data=lapop))

fit.party.ab.pooled <- svyglm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + party + bribe_victim*party, design = svydesign(ids=~1, strata=NULL, 
    weights=afro$wt, data=afro))
fit.party.ab.fe <- svyglm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + party + bribe_victim*party + as.factor(ccode_abc), 
    design = svydesign(ids=~1, strata=NULL, 
    weights=afro$wt, data=afro))

stargazer(fit.party.lapop.pooled, fit.party.lapop.fe, 
    fit.party.ab.pooled, fit.party.ab.fe,
  column.labels=c("LAPOP<br>Pooled", "LAPOP<br>FE", "Afrobarometer<br>Pooled", 
    "Afrobarometer<br>FE"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES", "Incumbent Copartisan", "Bribe*Copartisan"),
  notes = "Columns (1) through (4) include survey weights",
  keep.stat = c("n", "rsq"), omit="ccode_abc",
  type="html", out="pooledmodels_partisanship_bysurvey.html")

#############################
# APPENDIX H: VICTIMIZATION #
#############################

# APPENDIX TABLE 6

fit.p.victim.combined.pooled <- lm(bribe_victim ~ female + age + 
    edu + quantSES + party, data=combined.data.mlm)
fit.p.victim.combined.fe <- lm(bribe_victim ~ female + age + 
    edu + quantSES + party + as.factor(ccode_abc), data=combined.data.mlm)

stargazer(fit.p.victim.combined.pooled, fit.p.victim.combined.fe,
  se = starprep(fit.p.victim.combined.pooled, fit.p.victim.combined.fe),
  column.labels=c("Pooled", "Fixed Effects"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Female", "Age", "Education", 
    "SES", "Incumbent Copartisan"),
  keep.stat = c("n", "rsq"), omit="ccode_abc", notes="Robust standard errors",
  type="html", out="pooledmodels_bribevictim.html")

# APPENDIX TABLE 7

fit.party.control.combined.pooled <- lm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + party, data=combined.data.mlm)
fit.party.control.combined.fe <- lm(vote_incumbent ~ bribe_victim + female + age + 
    edu + quantSES + party + as.factor(ccode_abc), data=combined.data.mlm)

stargazer(fit.party.control.combined.pooled, fit.party.control.combined.fe,
  se = starprep(fit.party.control.combined.pooled, fit.party.control.combined.fe),
  column.labels=c("Pooled", "Fixed Effects"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES", "Incumbent Copartisan"),
  keep.stat = c("n", "rsq"), omit="ccode_abc", note="Robust standard errors",
  type="html", out="partisanshipcontrol.html")

########################################
# APPENDIX I: RESPONDENTS WITH CONTACT #
########################################

fit.ab.contact.pooled <- svyglm(paste("vote_incumbent ~", predictors.f), 
  design = svydesign(ids=~1, strata=NULL, weights=small.afro[small.afro$some_contact==1,]$wt, 
  data=small.afro[small.afro$some_contact==1,]))
fit.ab.contact.fe <- svyglm(paste("vote_incumbent ~", predictors.f, "+as.factor(ccode_abc)"), 
  design = svydesign(ids=~1, strata=NULL, weights=small.afro[small.afro$some_contact==1,]$wt, 
  data=small.afro[small.afro$some_contact==1,]))

stargazer(fit.ab.contact.pooled, fit.ab.contact.fe,
  column.labels=c("Pooled", "Fixed Effects"),
  dep.var.labels.include=FALSE,
  dep.var.caption="",
  model.names=FALSE,
  covariate.labels = c("Bribe Victim", "Female", "Age", "Education", 
    "SES"),
  keep.stat = c("n", "rsq"), omit="ccode_abc", notes="OLS regressions estimated with survey weights",
  type="html", out="ab_some_contact.html")

##############################################
#        COUNTRY-BY-COUNTRY RESULTS          #
# APPENDIX C: COUNTRY-LEVEL VOTING RESULTS   #
# APPENDIX E: COUNTRY-LEVEL TURNOUT RESULTS  #
# APPENDIX G: COUNTRY-LEVEL PARTISAN RESULTS #
##############################################

# Some Background Information on the Number of Countries in the Data
lapop.countries <- length(unique(lapop$ccode_abc))
afro.countries <- length(unique(afro$ccode_abc))
total.countries <- lapop.countries + afro.countries

# Create Collection Matrices
coefs.mat.voteinc <- matrix(NA, nrow=total.countries,
 ncol=2, dimnames=list(c(unique(lapop$ccode_abc), unique(afro$ccode_abc)), 
 c("beta.inc", "se.inc")))
coefs.mat.voteyn <- coefs.mat.voteinc
colnames(coefs.mat.voteyn) <- c("beta.yn", "se.yn")
coefs.mat.voteinc.party <- matrix(NA, nrow=total.countries,
 ncol=5, dimnames=list(c(unique(lapop$ccode_abc), unique(afro$ccode_abc)), 
 c("beta.bribe", "var.bribe", "beta.bribe.party", "var.bribe.party", "cov.bribe.party")))

# Run a loop that creates a subset of data by country, runs regressions with the two
 # different DVs of interest, and then collects the beta-hats and standard errors for 
 # the bribe victimization coefficient 
# We use the svydesign() and svyglm() commands from the survey package to estimate
 # models with survey weights

for (i in 1:total.countries){
 ifelse(i <= lapop.countries, 
  data <- svydesign(ids=~1, 
   strata=lapop[lapop$ccode_abc==unique(lapop$ccode_abc)[i],]$estratopri, 
   weights=lapop[lapop$ccode_abc==unique(lapop$ccode_abc)[i],]$wt,
   data=lapop[lapop$ccode_abc==unique(lapop$ccode_abc)[i],]),
  data <- svydesign(ids=~1,
   strata=NULL,
   weights=afro[afro$ccode_abc==unique(afro$ccode_abc)[i-lapop.countries],]$withinwt,
   data=afro[afro$ccode_abc==unique(afro$ccode_abc)[i-lapop.countries],]))
 
 fit.temp1 <- svyglm(vote_incumbent ~ bribe_victim + female + age + edu + quantSES, design = data)
 fit.temp2 <- svyglm(vote_yes ~ bribe_victim + female + age + edu + quantSES, design = data)
 fit.temp3 <- svyglm(vote_incumbent ~ bribe_victim + female + age + edu + quantSES + party +
  bribe_victim*party, design = data)

 coefs.mat.voteinc[i, "beta.inc"] <- coef(fit.temp1)["bribe_victim"]
 coefs.mat.voteinc[i, "se.inc"] <- sqrt(vcov(fit.temp1)["bribe_victim", "bribe_victim"])
 coefs.mat.voteyn[i, "beta.yn"] <- coef(fit.temp2)["bribe_victim"]
 coefs.mat.voteyn[i, "se.yn"] <- sqrt(vcov(fit.temp2)["bribe_victim", "bribe_victim"])
 coefs.mat.voteinc.party[i, "beta.bribe"] <- coef(fit.temp3)["bribe_victim"]
 coefs.mat.voteinc.party[i, "var.bribe"] <- vcov(fit.temp3)["bribe_victim","bribe_victim"]
 coefs.mat.voteinc.party[i, "beta.bribe.party"] <- coef(fit.temp3)["bribe_victim:party"]
 coefs.mat.voteinc.party[i, "var.bribe.party"] <- vcov(fit.temp3)["bribe_victim:party","bribe_victim:party"]
 coefs.mat.voteinc.party[i, "cov.bribe.party"] <- vcov(fit.temp3)["bribe_victim","bribe_victim:party"]
}

# APPENDIX FIGURE 1: Plot Coefficients by Magnitude

f1_res = data.frame(est = coefs.mat.voteinc[,"beta.inc"][order(coefs.mat.voteinc[,"beta.inc"])],
                 ses = coefs.mat.voteinc[,"se.inc"][order(coefs.mat.voteinc[,"beta.inc"])],
                 lab = rownames(coefs.mat.voteinc)[order(coefs.mat.voteinc[,"beta.inc"])])

f1_res %>%
  ggplot(aes(1:length(est),est)) +
  geom_point() +
  geom_errorbar(aes(ymin=est-2*ses,
                    ymax=est+2*ses),
                width=0) +
  geom_text(aes(label=lab),
            angle=90,
            nudge_y = .05 + 2*f1_res$ses,
            size=2.5) +
  geom_hline(yintercept=0,linetype=2) +
  scale_x_continuous(breaks=NULL) +
  labs(x="Country",
       y="Coefficient on Bribe Exposure",
       title="Effect of Bribe Exposure on Voting for Incumbent") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        plot.title = element_text(hjust=.5)) 

# Summary of Coefficients

coef.stats <- function(coef.mat){
 print("Total Number of Statistics in Matrix")
 number.est <- sum(is.na(coef.mat[,"beta.inc"])==0)
 print(number.est)
 print("Number of Negative Coefficients for Incumbent Voting")
 print(sum(coef.mat[,"beta.inc"]<0, na.rm=T))
 print("Proportion of Incumbent Voting Coefficients That are Negative")
 print(sum(coef.mat[,"beta.inc"]<0, na.rm=T) / number.est)
 print("Number of Negative Coefficients that are Significant at p < 0.05")
 print(sum(coef.mat[,"beta.inc"] < 0 & 
  2*pt(abs(coef.mat[,"beta.inc"])/coef.mat[,"se.inc"], df=Inf, lower=FALSE) 
  < 0.05, na.rm=T))
 print("Proportion of Negative Coefficients that are Significant at p < 0.05")
 print(sum(coef.mat[,"beta.inc"] < 0 & 
  2*pt(abs(coef.mat[,"beta.inc"])/coef.mat[,"se.inc"], df=Inf, lower=FALSE) 
  < 0.05, na.rm=T) / number.est)
 print("Number of Positive Coefficients that are Significant at p < 0.05")
 print(sum(coef.mat[,"beta.inc"] > 0 & 
  2*pt(abs(coef.mat[,"beta.inc"])/coef.mat[,"se.inc"], df=Inf, lower=FALSE) 
  < 0.05, na.rm=T))
 print("Proportion of Positive Coefficients that are Significant at p < 0.05")
 print(sum(coef.mat[,"beta.inc"] > 0 & 
  2*pt(abs(coef.mat[,"beta.inc"])/coef.mat[,"se.inc"], df=Inf, lower=FALSE) 
  < 0.05, na.rm=T) / number.est)
}

coef.stats(coefs.mat.voteinc)
coef.stats(coefs.mat.voteinc[1:lapop.countries,])
coef.stats(coefs.mat.voteinc[(lapop.countries+1):total.countries,])

median(coefs.mat.voteinc[,"beta.inc"])
median(coefs.mat.voteinc[1:lapop.countries,"beta.inc"])
median(coefs.mat.voteinc[(lapop.countries+1):total.countries,"beta.inc"])

# Table
coefs.mat.voteinc.2 <- round(cbind(coefs.mat.voteinc, 
 2*pt(abs(coefs.mat.voteinc[,"beta.inc"])/coefs.mat.voteinc[,"se.inc"], 
 df=Inf, lower=FALSE)), 3)
colnames(coefs.mat.voteinc.2) <- c("beta.inc", "se.inc", "p-val")
HTML(coefs.mat.voteinc.2, "coefs_country-by-country.html")

# APPENDIX FIGURE 2: Turnout

ad_res = data.frame(est = coefs.mat.voteyn[,"beta.yn"][order(coefs.mat.voteyn[,"beta.yn"])],
                 ses = coefs.mat.voteyn[,"se.yn"][order(coefs.mat.voteyn[,"beta.yn"])],
                 lab = rownames(coefs.mat.voteyn)[order(coefs.mat.voteyn[,"beta.yn"])])

ad_res %>%
  ggplot(aes(1:length(est),est)) +
  geom_point() +
  geom_errorbar(aes(ymin=est-2*ses,
                    ymax=est+2*ses),
                width=0) +
  geom_text(aes(label=lab),
            angle=90,
            nudge_y = .05 + 2*ad_res$ses,
            size=2.5) +
  geom_hline(yintercept=0,linetype=2) +
  scale_x_continuous(breaks=NULL) +
  labs(x="Country",
       y="Coefficient on Bribe Exposure",
       title="Effect of Bribe Exposure on Turnout") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        plot.title = element_text(hjust=.5)) 

# APPENDIX TABLE 4: Summary

# Putting Betas and SEs into Dataset with Country-Level Variables
country.combined <- merge(country, as.data.frame(coefs.mat.voteinc), by.x="ccode_abc", 
 by.y="row.names")

coef.stats(country.combined)
coef.stats(country.combined[country.combined$p_polity2 >= 6,])
coef.stats(country.combined[country.combined$p_polity2 < 6,])

# APPENDIX FIGURE 3: Plot by Polity Score

# Create Vectors for Plot
country.var <- country.combined$p_polity2
beta.inc.temp <- country.combined$beta.inc
se.inc.temp <- country.combined$se.inc
ccode_abc.temp <- country.combined$ccode_abc

# Plot
plot(country.var, beta.inc.temp, xlab="Polity Score", ylab="Coefficient on Bribe Exposure", 
 main="Effect of Bribe Exposure on Voting for the Incumbent\nBy Polity Score",
 type="n")
text(country.var, beta.inc.temp, labels=ccode_abc.temp)
segments(country.var, beta.inc.temp-1.96*se.inc.temp, country.var, beta.inc.temp+1.96*se.inc.temp)
abline(h=0, lty=2)
abline(lm(beta.inc.temp ~ country.var), lty=5)

# APPENDIX FIGURE 4: Partisanship

order.party <- order(coefs.mat.voteinc.party[,"beta.bribe"])

f2_res = data.frame(est = c(coefs.mat.voteinc.party[,"beta.bribe"][order.party],
                         coefs.mat.voteinc.party[,"beta.bribe"][order.party] + 
                         coefs.mat.voteinc.party[,"beta.bribe.party"][order.party]),
                 ses = c(sqrt(coefs.mat.voteinc.party[,"var.bribe"][order.party]),
                         sqrt(coefs.mat.voteinc.party[,"var.bribe"][order.party] +
                         coefs.mat.voteinc.party[,"var.bribe.party"][order.party] +
                         2*coefs.mat.voteinc.party[,"cov.bribe.party"][order.party])), 
                 typ = rep(c("Non-Copartisans","Copartisans"),
                           each=length(coefs.mat.voteinc.party[order.party])),
                 lab = c(paste(rownames(coefs.mat.voteinc.party)[order.party], "Non"),
                         paste(rownames(coefs.mat.voteinc.party)[order.party], "Co")))


ggplot(f2_res %>% filter(typ=="Non-Copartisans"),
       aes(1:length(est)-.25,est)) +
  geom_point(color='darkgrey') +
  geom_text(aes(label=lab),angle=90,
            nudge_y = .2 + 2*f2_res$ses[f2_res$typ=="Non-Copartisans"],
            color='darkgrey',
            size=2.5) + 
  geom_errorbar(aes(ymin=est-2*ses,
                    ymax=est+2*ses),
                width=0,
                color='darkgrey') +
  geom_point(data=f2_res %>% filter(typ=="Copartisans"),
             aes(1:length(est)+.25,est)) +
  geom_text(data=f2_res %>% filter(typ=="Copartisans"),
            aes(x=1:length(est)+.25,
                label=lab),angle=90,
            size=2.5,
            nudge_y = -.2 - 2*f2_res$ses[f2_res$typ=="Copartisans"]) +
  geom_errorbar(data=f2_res %>% filter(typ=="Copartisans"),
                aes(x = 1:length(est)+.25,
                    ymin=est-2*ses,
                    ymax=est+2*ses),
                width=0) +
  geom_hline(yintercept=0,linetype=2) +
  scale_x_continuous(breaks=NULL) +
  labs(x="Country (in Order of Magnitude of Non-Copartisan Coefficient)",
       y="Point Estimate",
       title="Effect of Bribe Exposure on Voting for Incumbent\nby Partisanship") +
  annotate("text",c(5,5),c(-.75,-.85),
           label=c("Non-Copartisans",
                   "Copartisans"),
           color=c("darkgrey","black")) +
  ylim(c(-1.5,.75)) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        plot.title = element_text(hjust=.5)) 

# Count of Significantly More Negative than Other Effect for Non-Copartisans
 # Interaction term is positive (i.e., copartisans > non-copartisans) and significant
sum(coefs.mat.voteinc.party[,"beta.bribe.party"] > 0 & 
 2*pnorm(-abs(coefs.mat.voteinc.party[,"beta.bribe.party"]/sqrt(coefs.mat.voteinc.party[,"var.bribe.party"])))<0.05)

# Count of Significantly More Positive than Other Effect for Non-Copartisans
 # Interaction term is negative (i.e., copartisans < non-copartisans) and significant
sum(coefs.mat.voteinc.party[,"beta.bribe.party"] < 0 & 
 2*pnorm(-abs(coefs.mat.voteinc.party[,"beta.bribe.party"]/sqrt(coefs.mat.voteinc.party[,"var.bribe.party"])))<0.05)


# End of File